function y = func_delta_from_nz(nz, nx, CONSTS)
    
    eta   = CONSTS.eta;
    eps_a = CONSTS.eps_a;

    nxy = func_nxy_from_nz(nz, CONSTS);
    ny  = func_ny_from_nxy(nz, nx, CONSTS);
    [n1 n2] = func_nk_from_nz(nz, CONSTS);
    [J1 J2 J1_1 J2_1] = func_J_k(nxy, nx, CONSTS);
    
    y = -(n2-n1).*((eps_a-nx.^2)./(eps_a*nxy.^2)) - ...
        n1.*((eta/eps_a)*J2.*J1_1+(1i*ny./(nxy.^2)).*(J2+(eta/eps_a)*J1_1)) + ...
        n2.*((eta/eps_a)*J1.*J2_1+(1i*ny./(nxy.^2)).*(J1+(eta/eps_a)*J2_1)) + ...
        nz.*(eta/(eps_a)).*(nx./(nxy.^2)).*(J1+J1_1-J2-J2_1);

end

